function[px,py]= period(x,y,ofac,fhi),

%   np = number of frequencies examined (= length of output freqs file)
%   ofac = oversampling parameter (typically >= 4 for best results)
%   hifac = input parameter defined as fhi/fc, where:
%   fhi is the highest frequency examined, and
%   fc is the Nyquist frequency (=N/(2T)).

%   nuhifac=questdlg('Enter fhi or hifac?','Cutoff Frequency','fhi','hifac','fhi');
%   if strcmp(nuhifac,'fhi')
% 	fhi = getfhi;
% 	hifac=fhi*2*(tmax-tmin)/n;
%   elseif strcmp(nuhifac,'hifac')
% 	hifac = gethifac;
% 	fhi=hifac*n/(2*(tmax-tmin));
% 	fprintf('\nfhi = %f (in Hertz)',fhi);
% 	info{lines}=sprintf('fhi = %f (in Hertz)',fhi);
% 	lines=lines+1;
%   end

tmin=min(x);
tmax=max(x);
n=length(x);

% ofac = getofac;
% pause(0.1);
% close(findobj('name','ofac'));
% 
% fprintf('\nhifac = %f',hifac);
% info{lines}=sprintf('hifac = %f',hifac);
% lines=lines+1;
% fprintf('\nofac = %f',ofac);
% info{lines}=sprintf('ofac = %f',ofac);
% lines=lines+1;

hifac=fhi*2*(tmax-tmin)/n;
np=ofac*hifac*n*0.5;
%nout=floor(0.5*ofac*hifac*n);
nout=round(0.5*ofac*hifac*n);
ave=mean(y);
variance=var(y);
xmin=tmin;
xmax=tmax;
xdif=xmax-xmin;
xave=0.5*(xmax+xmin);
pymax=0;
pnow=1/(xdif*ofac);
hbar = waitbar(0,'JVALS...');
for jval=1:n
   	waitbar(jval/n);

	arg=2*pi*((x(jval)-xave)*pnow);
	wpr(jval)=-2*sin(0.5*arg)^2;
	wpi(jval)=sin(arg);
	wr(jval)=cos(arg);
	wi(jval)=wpi(jval);
end;
close(hbar);
hbar = waitbar(0,'IVALS...');

for ival=1:nout
   	waitbar(ival/nout);
	px(ival)=pnow;
	sumsh=0;
	sumc=0;
	for jval=1:n
		c=wr(jval);
		s=wi(jval);
		sumsh=sumsh+s*c;
		sumc=sumc+(c-s)*(c+s);
	end
	wtau=0.5*atan2(2*sumsh,sumc);
	swtau=sin(wtau);
	cwtau=cos(wtau);
	sums=0;
	sumc=0;
	sumsy=0;
	sumcy=0;
	for jval=1:n
		s=wi(jval);
		c=wr(jval);
		ss=s*cwtau-c*swtau;
		cc=c*cwtau+s*swtau;
		sums=sums+ss^2;
		sumc=sumc+cc^2;
		yy=y(jval)-ave;
		sumsy=sumsy+yy*ss;
		sumcy=sumcy+yy*cc;
		wtemp=wr(jval);
		wr(jval)=(wr(jval)*wpr(jval)-wi(jval)*wpi(jval))+wr(jval);
		wi(jval)=(wi(jval)*wpr(jval)+wtemp*wpi(jval))+wi(jval);
	end
%	py(ival)=0.5*(sumcy^2/sumc+sumsy^2/sums)/variance;
		py(ival)=(sumcy^2/sumc+sumsy^2/sums)/n;  % changed 23.7.2004
	%WRITE OUTPUT
	freqs(ival,1)=px(ival);
	freqs(ival,2)=py(ival);
	pnow=pnow+1/(ofac*xdif);
end

close(hbar);
return